
%%% The Ambiguity of Fishing for Fun, October 2022
%%% This function evaluates the loglikelihood and its gradient for the CARA utility model

function[LnL , gradient]=likelihoodGrad_cara(K,Q,G,Z,P,Y,C,I,IJ,ID,parameters)

delta=[parameters(3), parameters(4)]'; %
vartheta=[parameters(8),parameters(9),parameters(10),parameters(11),parameters(12),parameters(13), parameters(14)];
u_W=zeros(length(Q),K);
Wk=zeros(length(Q),K); 

for k=1:K
   Wk(:,k)=parameters(1).*Q(:,k)+parameters(2).*G(:,k)+Z*delta; 
   u_W(:,k)=(1-exp(-parameters(5).*Wk(:,k)))/parameters(5); 
end

F=sum(u_W.*P,2)+parameters(6).*C; %Risk aversion for both kept and released fish
EF=exp(F);

%Partial derivatives with respect to theta and r (risk aversion) 
f_theta1=sum(((1-parameters(5)*u_W).*Q).*P,2);
f_theta2=sum(((1-parameters(5)*u_W).*G).*P,2);
f_theta3=Z(:,1).*sum((1-parameters(5)*u_W).*P,2);
f_theta4=Z(:,2).*sum((1-parameters(5)*u_W).*P,2);
derivativeAlpha=(((1/parameters(5))*(Wk+1/parameters(5))).*(1-parameters(5)*u_W))-1/(parameters(5)^2);
f_theta5=sum(derivativeAlpha.*P,2);

%Calculating the opt-out utility
Eopt_out=exp(parameters(7).*ones(length(C),1)+Y*vartheta'+parameters(6).*C);

%Separating trips A and B from the opt-out trip alternative
ef=EF.*(1-IJ);
eopt_out=Eopt_out.*IJ;
Denom=ef+eopt_out;
InactiveCosts=C.*IJ;
sumtheta1=f_theta1.*EF.*(1-IJ);
sumtheta2=f_theta2.*EF.*(1-IJ);
sumtheta3=f_theta3.*EF.*(1-IJ);
sumtheta4=f_theta4.*EF.*(1-IJ);
sumtheta5=f_theta5.*EF.*(1-IJ);
sumtheta6=C.*EF.*(1-IJ);

%Forming the denominator
sume=accumarray(ID, Denom); 
sumeInactiveCost=accumarray(ID, InactiveCosts);
sumActiveTrips=accumarray(ID, ef);
sumInactiveTrips=accumarray(ID,eopt_out);
Sumtheta1=accumarray(ID,sumtheta1);
Sumtheta2=accumarray(ID,sumtheta2);
Sumtheta3=accumarray(ID,sumtheta3);
Sumtheta4=accumarray(ID,sumtheta4);
Sumtheta5=accumarray(ID,sumtheta5);
Sumtheta6=accumarray(ID,sumtheta6);
SUME=kron(sume,[1 1 1]'); 
SUMEInactiveCost=kron(sumeInactiveCost,[1 1 1]');
SUMEActiveTrips=kron(sumActiveTrips,[1 1 1]');
SUMEInactiveTrips=kron(sumInactiveTrips,[1 1 1]');
SUMEtheta1=kron(Sumtheta1,[1 1 1]');
SUMEtheta2=kron(Sumtheta2,[1 1 1]');
SUMEtheta3=kron(Sumtheta3,[1 1 1]');
SUMEtheta4=kron(Sumtheta4,[1 1 1]');
SUMEtheta5=kron(Sumtheta5,[1 1 1]');
SUMEtheta6=kron(Sumtheta6,[1 1 1]');

F=ef./SUME; 
FTerm=F(F~=0);
FTermI=I(F~=0);
S=eopt_out./SUME; 
STerm=S(S~=0);
STermI=I(S~=0);

%Calculating lnL 
LnL=-(sum(log(FTerm).*FTermI)+sum(log(STerm).*STermI));

%Calculating the gradient of lnL
if nargout>1
    gradient(1)=-(sum(I.*(1-IJ).*(f_theta1-(SUMEtheta1./SUME)))-sum(I.*IJ.*(SUMEtheta1./SUME)));
    gradient(2)=-(sum(I.*(1-IJ).*(f_theta2-(SUMEtheta2./SUME)))-sum(I.*IJ.*(SUMEtheta2./SUME))); 
    gradient(3)=-(sum(I.*(1-IJ).*(f_theta3-(SUMEtheta3./SUME)))-sum(I.*IJ.*(SUMEtheta3./SUME))); 
    gradient(4)=-(sum(I.*(1-IJ).*(f_theta4-(SUMEtheta4./SUME)))-sum(I.*IJ.*(SUMEtheta4./SUME))); 
    gradient(5)=-(sum(I.*(1-IJ).*(f_theta5-(SUMEtheta5./SUME)))-sum(I.*IJ.*(SUMEtheta5./SUME)));  
    gradient(6)=-(sum(I.*(1-IJ).*(C-(SUMEtheta6+SUMEInactiveCost.*SUMEInactiveTrips)./SUME))+sum(I.*IJ.*(SUMEInactiveCost-(SUMEtheta6+SUMEInactiveCost.*SUMEInactiveTrips)./SUME))); 
    gradient(7)=-(sum(I.*IJ.*(SUMEActiveTrips./SUME))-sum(I.*(1-IJ).*(SUMEInactiveTrips./SUME))); 
    gradient(8)=-(sum(I.*IJ.*Y(:,1).*(SUMEActiveTrips./SUME))-sum(I.*(1-IJ).*Y(:,1).*(SUMEInactiveTrips./SUME))); 
    gradient(9)=-(sum(I.*IJ.*Y(:,2).*(SUMEActiveTrips./SUME))-sum(I.*(1-IJ).*Y(:,2).*(SUMEInactiveTrips./SUME))); 
    gradient(10)=-(sum(I.*IJ.*Y(:,3).*(SUMEActiveTrips./SUME))-sum(I.*(1-IJ).*Y(:,3).*(SUMEInactiveTrips./SUME))); 
    gradient(11)=-(sum(I.*IJ.*Y(:,4).*(SUMEActiveTrips./SUME))-sum(I.*(1-IJ).*Y(:,4).*(SUMEInactiveTrips./SUME))); 
    gradient(12)=-(sum(I.*IJ.*Y(:,5).*(SUMEActiveTrips./SUME))-sum(I.*(1-IJ).*Y(:,5).*(SUMEInactiveTrips./SUME))); 
    gradient(13)=-(sum(I.*IJ.*Y(:,6).*(SUMEActiveTrips./SUME))-sum(I.*(1-IJ).*Y(:,6).*(SUMEInactiveTrips./SUME))); 
    gradient(14)=-(sum(I.*IJ.*Y(:,7).*(SUMEActiveTrips./SUME))-sum(I.*(1-IJ).*Y(:,7).*(SUMEInactiveTrips./SUME)));
end
